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1 Introduction 

Wavelets are a powerful new mathematical tool which offers the possibility to treat in a 
natural way quantities characterized by several length scales. In this article we will show 
how wavelets can be used to solve partial differential equations which are characterized by 
widely varying length scales and which are therefore hardly accessible by other numerical 
methods. The standard way to solve partial differential equations is to express the solution 
as a linear combination of so-called basis functions. These basis functions can for instance 
be plane waves, Gaussians or finite elements. Having discretized the differential equation 
in this way makes it amenable to a numerical solution. Wavelets are just another basis 
set which however offers considerable advantages over alternative basis sets. Its main 
advantages are: 

1. The basis set can be improved in a systematic way: 

If one wants the solution of the differential equation with higher accuracy one can 
just add more wavelets in the expansion of the solution. This will not lead to any 
numerical instabilities. 

2. Different resolutions can be used in different regions of space: 

If the solution of the differential equation is varying particularly rapidly in a par- 
ticular region of space one can increase the resolution in this region by adding more 
high resolution wavelets centered around this region. 

3. There are few topological constraints for increased resolution regions: 

The regions of increased resolution can be chosen in arbitrarily, the only requirement 
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being that a region of higher resolution be contained in a region of the next lower 
resolution. 

4. The matrix elements of the differential operators are very easy to calculate 

5. The numerical effort scales linearly with respect to system size: 
Three-dimensional problems of realistic size require usually a very large number 
of basis functions. It is therefore of utmost importance, that the numerical effort 
scales only linearly (and not quadratically or cubically) with respect to the number 
of basis functions. If one uses iterative matrix techniques, this requirement is equiv- 
alent to the two requirements, namely that the matrix vector multiplications which 
are necessary for all iterative methods can be done with linear scaling and that the 
number of matrix vector multiplications is independent of the problem size. The 
first requirement is fulfilled since the matrix representing the differential operator 
is sparse. The second requirement is related to the availability of a good precon- 
ditioning scheme which can be easily found by analyzing the Fourier properties of 
wavelets. 

2 A first tour of some wavelet families 

Many families of wavelets have been proposed in the mathematical literature. If one wants 
to use wavelets for the solution of differential equations, one therefore has to choose one 
specific family which is most advantageous for the intended application. Within one family 
there are also members of different degree. We believe that the so-called bi-orthogonal 
interpolating wavelets H are the most useful ones in the context of differential equations 
and we will therefore mainly concentrate on this class. Each wavelet family is characterized 
by two functions, the mother scaling function <ft and the mother wavelet ip. For the case 
of a fourth order interpolating wavelet they are shown in Figure |l|. 
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Figure 1: The interpolating scaling function and wavelet of degree 4 

Another family which will be introduced is the Haar wavelet family shown in Figure 
It is too crude to be useful for any numerical work, but its simplicity will help us to 
illustrate some basic wavelet concepts. 
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Figure 2: The Haar scaling function and wavelet. 



To obtain a basis set at a certain resolution level k one can use all the integer trans- 
lations of the mother scaling function of some wavelet family. 

<f> k i {x) = <f ) {2 k x-i) (1) 

Note that with this convention higher resolution corresponds to larger values of k. 
Exactly the same scaling and shifting operations can of course also be applied to the 
wavelets. 

^{x)=^{2 k x-i) (2) 

This set of wavelet basis functions can be added as a basis to the scaling functions as will 
be explained in the following for the case of the Haar wavelet family. 

3 The Haar wavelet 

In the case of the Haar family, any function which can exactly be represented at any level 
of resolution is necessarily piecewise constant. One such function is shown in Figure |3|. 
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Figure 3: function at resolution level 4. 



Evidently this function can be written as a linear combination of the scaling functions 
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/ = E*i#(*) ( 3 ) 

i=0 
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where sf = /(i/16). 

Another, more interesting, possibility consists of expanding a function with respect to 
wavelets of different resolution. This is possible because a scaling function (and wavelet) 
at resolution level k is always a linear combination of a scaling function and a wavelet at 
the next coarser level k — 1 as shown in Figure 
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Figure 4: A skinny scaling function is a linear combination of a fat scaling function and 
a wavelet. 



Using this relation, we can write any linear combination of the two scaling functions 
4>2i(x) and 4>2i+i( x ) as a linear combination of and ?/>f _1 (a;). 

Denoting the expansion coefficients with respect to ipf(x) as d\, we obviously obtain 

„fc-l _ I A , I ,* . jfc-1 _ Ifc _ I fc / 4 \ 

So to calculate the expansion coefficients with respect to the scaling functions at the 
next coarser level, we have to take an average over expansion coefficients at the higher 
resolution level. Because we have to take some weighted sum these coefficients are denoted 
by s. To get the expansion coefficients with respect to the wavelet, we have to take some 
weighted difference and the coefficients are accordingly denoted by d. The wavelet part 
contains mainly high frequency components and by doing this transformation we therefore 
peel off the highly oscillatory parts of the function. The remaining part represented by 
the coefficients ~ 1 is therefore smoother. For the case of our example in Figure |3] the 
remaining scaling function part after one transformation step is shown in Figure [5]. 

For any data set whose size is a power of 2, we can now apply this transformation 
repeatedly. In each step the number of s coefficients will be cut into half. So we have to 
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Figure 5: The function from Figure |3| at resolution level 3. 



stop the procedure as soon as there is only one s coefficient left. Such a series of trans- 
formation steps is called a forward Haar wavelet transform. The wavelet representation 
of the function in Equation [| is then 

/ = + dSVtfOn) + £<W(*) • ( 5 ) 

i=l i=l 8=1 

Note that in both cases we need exactly 16 coefficients to represent the function. 
Functional representations of this type will be the focus of this article. 

By doing a backward wavelet transform, we can go back to the original expansion 
of Equation [| Starting at the lowest resolution level, we have to split up each scaling 
function and wavelet on the coarse level into scaling functions at the finer level. 

= + ; 4ft = ( 6 ) 
4 The concept of Mult i- Resolution Analysis 

In the previous sections a very intuitive introduction to wavelet theory was given. The 
formal theory behind wavelets is called Multi- Resolution Analysis (MRA). The reader 
interested in the formal theory can consult Daubechies book. We will list here only a few 
facts which are useful for numerical work. 

A bi-orthogonal wavelet family of degree m is characterized by 4 finite filters denoted 
by hj, hj, gj, gj. A filter is just a short vector which is used in convolutions. Those filters 
satisfy certain orthogonality and symmetry relations. Scaling functions and wavelets at 
a coarse level can be written as the linear combinations of scaling functions at a higher 
resolution level. These important relations are called refinement relations. 

m 

<j>(x) = £ h j( f>(2x-j) (7) 

j=-m 
m 

^0*0 = E 9j<P(2x-j) (8) 

j=-m 
m 

4>{x) = 2 £ M(2x-j) (9) 



${x) = 2 £ gj^x-j) (10) 

j=-m 

The expansion coefficients at different resolution levels are related by the wavelet 
transform equations. The analysis (forward) transform is given by 

m m 



E ; dt 1 = E (ii) 



j=—m j=—m 

and a wavelet synthesis (backward) transform is given by 

m/2 m/2 

4 +1 = E h 2j stj + g 2j dtj ; 4+\= £ /,,,-. : .s? j • Uu^'l'i , (12) 

j=-m/2 j=-m/2 

These two equations are generalizations of Equations (|j) and (|6|) which we derived in an 
intuitive way and with a different normalization convention. 

The fundamental functions satisfy the following orthogonality relations 

^{x)<f) k j {x)dx = 5 id (13) 

$(x)(j)](x)dx = 0, k > q (14) 

^(x)4>](x)dx = 0, k > q (15) 

ipf(x)'ip](x)dx = 5 k:g S i:j (16) 

The scaling function is usually normalized to 1 

{x)dx = 1 (17) 

5 The fast wavelet transform 

Let us first look at the forward transform given by Equation [TTJ . The peeling off of the 
high frequency components in the forward transform can be illustrated in the following 

way: 



s 4 — > s 3 — > s 2 — > s 1 — > s° 



d 3 d 2 d 1 d° 



We note that just two arrays of length n (where n is a power of 2) are necessary to do 
the transform as shown below: 
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original data 

Array 1: s 4 , sf s\ sj s\ s\ 

after first sweep 

Array 2: s[j s i s l s l s l s l 

after second sweep 
Array 1: Sg S2 s| c?q 

after third sweep 

Array 2: sj cfj d 2 Q d\ 

final data 

Array 1: s° d Q Q d l d\ d 2 Q d\ 



4444444444 

b Q b 7 b 8 b g b 1Q b n b 12 b 13 b 1A b 15 

Sq S'j d/Q d-^ d<2 d^ d^ d^ d^ d^ 

^"2 ^3 ^0 d-^ d<2 d^ d^ d^ d^ 

d<2 d^ ^0 d^ d<2 dr^ d^ d^ d^ 

^"2 ^3 ^0 d-^ d<2 d^ d^ d?j d'j 



Note that this transformation from the "original data" to the "final data" corresponds 
exactly to the transformation done in an intuitive way to get from Equation [5] to Equa- 
tion |5|. Just as in the case of a Fast Fourier transform we have Log2(n) sweeps to do 
a full transform. However in the case of the wavelet transform the active data set (the 
s coefficients) is cut into half in each sweep. If our filters h and g have length 2m the 
operation count is then given by 2m(n + n/2 + n/4 + ...). Replacing the finite geometric 
series by its infinite value, the total operation count is thus given by 4mn 

The backward transform (Equation 0) can pictorially be represented by the following 
diagram: 



s 4 <— s 3 <— b 2 <— s 1 <— b° 



d 3 d 2 d 1 d° 



As can easily been seen the operation count is again 4mn and again it can be done 
with 2 arrays of length n. Since each sweep in a wavelet transform is a linear operation it 
can be represented by a matrix. Denoting the matrix for one sweep in a forward transform 
by F and in a backward transform by B we have 

F T = F- 1 = B ; B T = B- 1 = F (18) 

where the tilde on the matrix means that the filter coefficients necessary to fill the matrix 
are replaced by their dual counterparts. Obviously all these matrices are sparse and 
banded. 

Backward wavelet transforms can also be used to make plots of scaling functions and 
wavelets. To generate the scaling function we start with a data set where Sq = 1 and 
d\ = for all possible i's and /c's up to a maximum resolution level k = K. In the wavelet 
case the initial data set is Sq = 0, d® = 1, and d\ = for all other values of i and k up 
to the maximal resolution K. By doing repeated backward transform sweeps, we express 
these two functions by skinnier and skinnier scaling functions and the s coefficients will 
finally be the functional values within the resolution of the eye. 
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6 Interpolating wavelets 



In addition to being advantageous as basis sets, interpolating wavelets are also conceptu- 
ally the simplest wavelets and we will therefore briefly describe their construction. The 
construction of interpolating wavelets is closely connected to the question of how to con- 
struct a continuous function f(x) if only its values fa on a finite number of grid points i are 
known. One way to do this is by recursive interpolation. In a first step we interpolate the 
functional values on all the midpoints by using for instance the values of two grid points 
to the right and of two grid points to the left of the midpoint. These four functional 
values actually allow us to construct a third order polynomial and we can then evaluate 
it at the midpoint. In the next step, we take this new data set, which is now twice as 
large as the original one, as the input for a new midpoint interpolation procedure. This 
can be done recursively ad infinitum until we have a quasi continuous function. 

Let us now show, how this interpolation prescription leads to a set of basis functions. 
Denoting by the Kronecker a data set which has a nonzero entry only at the j-th 
position, we can write any initial data set also as a linear combination of such Kronecker 
data sets: fa = J2j fj$i-j- Now the whole interpolation procedure is clearly linear, i.e. the 
sum of two interpolated values of two functions is equal to the interpolated value of the 
sum of these two functions. This means that we can instead also take all the Kronecker 
data sets as the input for separate interpolation procedures, to obtain a set of functions 
4>(x — j). The final interpolated function is then identical to 

/(•'') Y.IM-'--./) (19) 

3 

If the initial grid values fa were the functional values of a polynomial of degree less 
than four, we obviously will have exactly reconstructed the original function from its 
values on the grid points. Since any smooth function can locally be well approximated 
by a polynomial, these functions <p(x) are good basis functions also in the case where / is 
not a polynomial and we will use them as scaling functions to construct a wavelet family. 

The first construction steps of an interpolating scaling function are shown below for 
the case of linear interpolation. The initial Kronecker data set is denoted by the big 
dots. The additional data points obtained after the first interpolation step are denoted 
by medium size dots and the additional data points obtained after the second step by 
small dots. 
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Continuing this process ad infinitum will then result in the function shown in the left 
panel of Figure |6|. If an higher order interpolation scheme is used the function shown in 
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the right panel of Figure |5] is obtained. 
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Figure 6: A Kronecker delta interpolated ad infinitum with linear interpolation (left 
panel) an 7-th order interpolation (right panel) . 



By construction it is clear, that <f)(x) has compact support. If an (m — l)-th order 
interpolation scheme is used, the filter length is (m — 1) and the support interval of the 
scaling function is [— (m — 1); (m — 1)]. 

It is also not difficult to see that the functions 4>{x) satisfy the refinement relation. 
Let us again consider the interpolation ad infinitum of a Kronecker data set which has 
everywhere zero entries except at the origin. We can now split up this process into 
the first step, where we calculate the half-integer grid point values, and a remaining 
series of separate ad infinitum interpolations for all half-integer Kronecker data sets, 
which are necessary to represent the data set obtained by the first step. Doing the 
ad-infinitum interpolation for a half integer Kronecker data set with a unit entry at 
position j, we obviously obtain the same scaling function, just compressed by a factor of 
2, (p(2x — j). If we are using a (m — l)-th order interpolation scheme (i.e. m input data 
for the interpolation process) we thus get the relation 

TO— 1 

4(x)= E 0(j/2)0(2x-j) (20) 

j=— m+l 

Comparing this equation with the refinement relation Equation [j] we can identify the 
first filter h as 

hj = 0(j/2) , j = -m + 1, m - 1 

For the case of third order interpolation the numerical values of h follow from the 
standard interpolation formula and are given by { -1/16 , , 9/16 , 1 , 9/16 , , -1/16 }. 

Let us next determine the filter h. Let us consider a function f(x) which is band- 
limited in the wavelet sense, i.e which can exactly be represented by a superposition of 
scaling functions at a certain resolution level K. 

/(*) = 5>f0f(*) 
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It then follows from the orthogonality relation Equation O that 

- f 4>f(x) f(x)dx (21) 
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Now we have seen above that with respect to interpolating scaling functions, a band- 
limited function is just any polynomial of degree less than or equal to m — 1, and that 
in this case the expansion coefficients are just the functional values at the grid points 
(Equation [T9]). We therefore have 



sf = / 4>f(x)f(x)dx = f 3 (22) 



which shows that the dual scaling function is the delta function. 

4>(x) = 5(x) (23) 

Obviously the delta function satisfies a trivial refinement relation S(x) = 25(2x) and from 
Equation |9] we conclude that hj = Sj From the symmetry relations for the filters the two 
remaining filters g(i) and g(i) can be determined and we have thus completely specified 
our wavelet family. 

Using these filters we can then determine the wavelet ip and its dual counterpart ip 
which turn out to be 

i){x) = <p{2x - 1) (24) 

(25) 

We see that the interpolating wavelet is a very special case in that its scaling function 
and wavelet have the same functional form and that the dual functions are related to the 
delta function. The non-dual functions are shown in Figure |I[ 

Lifting || is a very useful technique to modify an existing family of wavelets to meet 
specific needs. We can for instance lift the interpolating wavelets to obtain a new family 
whose wavelet has more vanishing moments M\. 



Mi — J ipf(x) x l dx 
which will for instance improve the frequency properties of the wavelet. 



7 Expanding functions in a wavelet basis 

As was demonstrated in the case of the Haar wavelet, there are two possible representa- 
tions of a function within the framework of wavelet theory. The first one is called scaling 
function representation and involves only scaling functions. The second is called wavelet 
representation and involves wavelets as well as scaling functions. Both representations 
are completely equivalent and exactly the same number of coefficients are needed in the 
case where one has uniform resolution. 
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The scaling function representation is given by 



/(*) = E «. 



Kmax iKmax , 



X 



(26) 



The coefficients Sj max can be calculated by integration through Equation [22|. Once 
we have a set of coefficients sf max we can use a full forward wavelet transform to obtain 
the coefficients of the wavelet representation 



Kmin iKmin, 



X) 



Kmax 

+ E E^fW 

j K=Kmin j 

Alternatively, one could also directly calculate the d coefficients by integration 

df = f $f{x) f(x)dx 



(27) 



(2f 



Equation 123] follows from the orthogonality relations [14] to ITB 



So we see that if we want to expand a function either in scaling functions or wavelets, 
we have to perform integrations at some point to calculate the coefficients. For general 
wavelet families this integration can be fairly cumbersome || and require especially in 2 
and 3 dimensions a substantial number of integration points. Furthermore it is not obvious 
how to do the integration if the function is only given in tabulated form. The interpolating 
wavelets discussed above are the glorious exception. Since the dual scaling function is a 
delta function ([23] ) and since the dual wavelet is a sum of delta functions (|25] ), one or a 
few data points are sufficient to do the integration exactly. One will therefore get exactly 
the same number of coefficients as one has data points and one has an invertible one-to-one 
mapping between the functional values on the grid and the expansion coefficients. This 
is even true in the case of nonuniform data sets, where we necessarily have to calculate 
the s and d coefficients directly by integration using [28j. As follows from Equation |23 



and |25|, one just needs the functional values at the data point at which the wavelet will 
be centered and a few data points at one lower resolution level around this center. If 
one wants to calculate the interpolating wavelet center at the high resolution grid point 
indicated by the fat arrow in the figure below, one needs in the case of the 4-th order 
interpolating wavelets the 4 additional points indicated by thin arrows which belong to a 
more coarse grid and are therefore always available even if the fine grid does not extend 
into this region. 



In the case where one wants to represent functions with several length scales which 
need inhomogeneous real space grid structure the wavelet representation allows a much 
more compact representation than the scaling function representation, since on can neglect 
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all the tiny d coefficients in the regions where one has little variation. To illustrate this 
let us look at the function / 

8 

f( x ) = J2 ex p(-( xl ) 2 ) 
1=1 

Evidently this function exhibits 8 different length scales. If one expands one simple Gaus- 
sian exp(— x 2 ) with respect to 4-th order interpolating scaling functions with a resolution 
of 1/16, one gets a reasonably small error of 10~ 6 . For the multi-scale function /, this 
error increases to more than 10~ 2 with the same resolution. If one however uses a scheme 
where one uses 32 wavelets on additional 5 resolution levels to improve the resolution as 
one approaches the origin one can again represent the function with an error of roughly 
10~ 6 ( it turns out that the expected 8 additional levels are not all needed). The total 
number of coefficients needed to represent the function in the interval [—2; 2] is then 4 x 
16 coefficients for the equal resolution (1/16) scaling function part plus 5 x 32 coefficients 
for the resolution enhancement with the wavelets, which makes all together 224 coeffi- 
cients. This has to be compared with the 1024 scaling function coefficients which would 
be needed to represent the function over the whole interval with the maximum resolution 
of (1/256), which we have obtained around the origin with this data compression scheme. 

8 Wavelets in 2 and 3 dimensions 

The easiest way to construct a wavelet basis in higher dimensional spaces is by form- 
ing product functions ||. For simplicity of notation we will only consider here the 2- 
dimensional case, the generalization to higher dimensional spaces being obvious. 
The space of all scaling functions of resolution level k is given by 

<eO»M/) = A (29) 
The wavelets consist of three types of products 

1>[8d& ta (x,y)=<ft l {z)ilb(y) (30) 

i>[ds)iMx,y)=i>U x )<f>i2(y) (3i) 

1>[dd\t Lja (x,y) = 1&(x)l4 l (y) (32) 

A wavelet transform step in the 2-dimensional setting is done by first transforming 
along the x and then along the y direction (or vice versa). 

9 The standard operator form 

In a bi-orthogonal wavelet basis it is natural to solve a differential equation in the colloca- 
tion sense. Let us recall that in the collocation method one has two functional spaces, the 
space of the basis function which are used to represent the solution and the space of the 
test functions which are used to multiply the differential equation from the left to obtain 
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a linear system of equations. In our case the expansion set are the scaling functions and 
wavelets while the test set are their dual counterparts. Lets consider the case of Poisson's 
equation 

VV = -4vrp. (33) 
Given the expansion of the charge density p in a wavelet basis 

Kmax 

p(*) = £sf mm 0f mm (*) + £ £cfVf(*) (34) 

j K=Kmin j 

we are looking for the wavelet expansion coefficients of the potential V. 

Kmax 

V(x) = J2Sf rnAn 4>f min (x)+ £ Y i Df^f{x) (35) 

j K=Kmin j 

Plugging in the expansion for p and V ( |4j) and ( |35|) in Poissons equation [33] and 
multiplying from the left with the dual wavelet collocation test space we obtain a system 
of equations 

A s v = p (36) 

where v is the vector containing both the s and d coefficients of the potential and p is 
the same vector for the charge density p. The matrix A s represents the Laplacian in this 
wavelet basis and one says that it has standard form. This standard form is graphically 
shown in Figure [7] . 

S(0) D(0) D(-l) D(-2) 

S(0) 



D(0) 



D(-l) 



D(-2) 



Figure 7: The structure of a matrix in the standard form. 

The problem with the standard form is that it is first of all rather complicated. There 
is coupling between all resolution levels and one has to calculate many different types of 
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matrix elements corresponding to all possible products of wavelets and scaling functions 
at different resolution levels and positions. The second point is that there are many blocks 
in that matrix which have no or only few zeroes. Let us look at the blocks representing 
the coupling between the scaling functions at the highest resolution level and the wavelets 
at the different resolution. In general each scaling function will extend over the whole 
computational volume and will therefore overlap with all the wavelets at any position. All 
these blocks will consequently have nonzero entries only. So this standard matrix form 
has more nonzero entries than we would like to have for optimal efficiency in the matrix 
vector multiplications which are required for all iterative linear equation solvers. 



10 The non-standard operator form 

The so-called nonstandard || form gives a much easier and efficient representation of our 
matrix. To derive it let us first assume, that our potential V and charge p are given in a 
scaling function basis. The Laplacian is then represented by a matrix A whose elements 
Aij are given by J (x)V 2 0^ (x)dx 
The matrix equation 

Av = p 

can graphically be represented in the following way: 



<sus> 

\ 

\ 




s 




s 



Now we can of course perform one step of a forward wavelet transform on all our 
data, i.e. both on the vector to be multiplied with the matrix and on the vector which 
is the result of this matrix times vector multiplication. Correspondingly we have then to 
transform the matrix A using the matrices whose properties are given in Equation [Tj| 

Fp = (FAF T )(Fv) 

Graphically this gives: 



s 

S 

<Slt§> 


\ 

\ 

<Slfl}> 


<DN§> 


:DHD> 



If we recursively applied wavelet transform to the upper S part we would obviously 
obtain the standard operator form. To get the nonstandard form, we have to add another 
step where we artificially enlarge our matrix A by putting in 5 blocks of zeroes as shown 
below: 
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<SIS> 














<SID> 





<DIS> 


<DH5> 



We see that our input and output vectors v and p also have to be adapted to this 
matrix structure leading to a redundant copy of the S data set. 

We can now recursively apply this 2-step procedure on the < S\S > block of the result- 
ing matrices. Doing this we obtain the so called non-standard form, which is graphically 
visualized in Figure || 





Figure 8: The structure of a matrix in the nonstandard form. 



As we see, we have now completely decoupled different resolution levels, since there 
are no blocks in this matrix between different levels. The coupling between different 
levels just enters trough the wavelet transforms which have to be interleaved with the 
application of this nonstandard operator form. We also see that all the nonzero blocks 
of this nonstandard matrix representation are strictly banded and the application of this 
matrix to a vector scales therefore linearly. 

The structure of the matrix in Figure |] is primarily valid for the case of uniform 
resolution where all the possible d coefficients at the highest resolution level are nonzero. 
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It can however easily be seen that this nonstandard form retains its advantage in a case of 
varying resolution where only some of the d coefficients are nonzero. If the nonredundant 
input data set is sparse, the redundant input data set will be sparse as well. Since all 
the blocks are banded, the redundant output set will be sparse as well. Finally the 
nonredundant output set will then be sparse as well. 

11 Calculation of differential operators in a wavelet 
basis 

As we have seen in the preceeding chapter we need the matrix elements 



(37) 
(38) 
(39) 
(40) 

for the application of an operator in the nonstandard form. Matrix elements on differ- 
ent resolution levels are related by simple scaling relations. So we just have to calculate 
these 4 matrix elements for one resolution level. On a certain resolution level, we can 
use the refinement relations to express the matrix elements involving wavelets in terms of 
matrix elements involving scaling functions (at a better resolution level) only. So we just 
have to calculate the basic integral aj 

Ql 

(x)—r(j)(x — i)dx (41) 
ox 1 

Using the refinement relations Equations |7] and for and we obtain 

Ql 

(x)——r(j)(x — i)dx (42) 
ox 1 



jut- 






dx 


/« 






dx 


/#( 


Ql 




dx 








dx 



y ' 2h u h„ I <p(2x — u)——f(p(2x — 2i — u)dx 

J2 Kh„2 l / (j){y)—(j){y -2i-fM+ v)dy 

(43) 
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We thus have to find the eigenvector a associated with the eigenvalue of 2 . 

= Q) <k (44) 

where the matrix A it j is given by 

Aid = hvh^ bj^i-v+p, (45) 

As it stands this eigensystem has a solution only if the rang of the matrix A — 2~ l I 
is less than its dimension. For a well defined differential operator, i.e if I is less than the 
degree of smoothness of the scaling function this will be the case. 



The system of equations |44] determines the a/s only up to a normalization factor. 
For the case of interpolating wavelets the normalization condition is easily found from the 
requirement that one obtains the correct result for the function x l . From the normalization 
of the scaling function (O) and from elementary calculus, it follows that 

<t>(x)^-;iJdx = I <j>(x)l\dx = l\ (46) 
ox 1 J 

On the other hand we know, that we can expand any polynomial of low enough degree 
exactly with the interpolating polynomials. The expansion coefficients are just i l by 



Equation ^2[ So we obtain 

d< 



^E'^-i) = £* , a i (47) 



By comparing Equation [46] and [47] we thus obtain the normalization condition 

y £i l a i = l\ (48) 

i 

The interpolating wavelet family offers also an important advantage for the calcu- 
lation of differential operators. Whereas in general derivative filters extend over the 
interval [—2m; 2m] their effective filter length is only [— m + 2; m — 2]. Since higher- 
dimensional wavelets are products of one- dimensional ones differential operators in the 
higher- dimensional case can easily be derived from the one-dimensional results. 

The standard operator form can not only be used for the application of differential 
operators, but also for other operations. If one want to transform for instance from one 
wavelet family (f> to another wavelet family $ the basic integral becomes 

Of = / $(x)(p(x - i)dx = Huh^i-u+fj, (49) 

Another use is for scalar products where the fundamental integral is 

ai = J (f)(x)4>(x - i)dx = ^ Kh^i-u+n (50) 
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12 Solving Poisson's equation for the U2 dimer 



Poisson's equation is a prototype differential equation and we want to solve it therefore 
as an illustration of wavelet theory. To demonstrate the power of the wavelet method we 
applied it to the most difficult system we could think of in the area of electronic structure 
calculations, namely the calculation of the electrostatic potential of a three dimensional 
U2 dimer [ IDfl . In this example, we clearly find widely varying length scales. The valence 
electrons have an extension of 5 atomic units, the Is core electrons of 2/100 atomic units 
and the nucleus itself was represented by a charge distribution with an extension of 1/2000 
atomic units. So all together the length scales varied by 4 orders of magnitude and two 
regions of increasing resolution (around each nucleus) were needed. In order to have quasi 
perfect natural boundary conditions we embedded the molecule in a computational volume 
of side length 10 4 atomic units. All together this necessitated 22 levels of resolution. Even 
though the potential itself varies by many orders of magnitude, we were able to calculate 
the solution with typically 7 digits of accuracy. We believe that it would not be possible 
with any other method to solve this kind of benchmark problem. 

The solution of Poisson's equation consists of several steps. Initially we have to find 
the wavelet expansion for a data set on a nonuniform real space grid structure shown in 
Figure |9] which represents the charge density. The resolution needed can in this exam- 
ple be estimated from the known extension and variation of the different atomic shells. 
Analogously to the one-dimensional case, this expansion can also easily be obtained for 
higher dimensional interpolating wavelets since all the dual function are related to delta 
functions. Let us point out, that also in this case the mapping from real space represen- 
tation to the wavelet representation is invertible, and we could thus get back exactly the 
same real space values if we evaluate the wavelet expansion on the grid points. 



Figure 9: A grid with two centers of increased resolution around the two nuclei. Only 3 
of the 22 levels used in the calculation are shown in this projection on a plane. 



Next we start a iteration loop for the potential. First we have to apply the Laplace 
operator to an approximate potential using the non-standard operator form. Subtract- 
ing from this result the charge density gives the residue vector which is the basis for all 
iterative methods @, such as steepest descent and conjugate gradient methods. Unfor- 
tunately the condition number of the Laplace matrix worsens when more high resolution 
levels are added and the number of iterations needed to obtain convergence would dramat- 
ically increase if we used straightforward iterative methods. It this therefore absolutely 
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necessary to use a preconditioned iterative method which will give a condition number 
which is independent of the maximal resolution. In a preconditioning scheme one has 
to find an approximate inverse matrix of the Laplace matrix. If the Laplace matrix is 
strongly diagonally dominant, then just the inverse of the diagonal part (which is again 
diagonal) will be a good approximate inverse. Whether the Laplace matrix is strongly 
diagonally dominant depends on the kind of wavelet family which is used. In a plane wave 
representation the Laplace matrix is strictly diagonal. If therefore our wavelet family has 
good frequency localization properties the resulting matrix will be strongly diagonally 
dominant. Unfortunately our favorite interpolating wavelets have a very poor frequency 
localization making an iterative solution practically impossible. It is therefore necessary to 
do the preconditioning step within another family such as the lifted interpolating wavelets 
which have much better frequency localization properties as shown in Figure |10[ Their 
improved frequence localization properties is related to the fact that several moments of 
the wavelet vanish. 



resolution levels. One has reasonable frequency separation in the lifted but not in the 
unlifted case. 

As discussed above the transformation into another wavelet family can also be done 
with the help of the non-standard operator form. The preconditioned residue vector is then 
used to update the potential and we go back to the beginning of the iteration. Using lifted 
interpolating wavelets with 2 vanishing moments we were able to reduce the norm of the 
residue vector by one order of magnitude with 3 iterations independent of the maximal 
resolution. Despite their poor frequency localization properties, unlifted interpolating 
wavelets have recently also been proposed for the solution of Poisson's equation ||. 

13 Outlook and conclusions 

Since we used mainly interpolating wavelets, all we did was essentially interpolating, 
which is one of the oldest technique in numerical analysis. However the framework pro- 




\ skinny 
fat 
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vided by wavelet theory puts this whole interpolation procedure on the new and powerful 
basis of multi-resolution analysis, expanding thus considerably the scope of interpolation 
based techniques. In particular it assigns basis functions to certain interpolation schemes. 
Wavelet based techniques allow us thus to solve differential equations which have several 
length scales and to do this with linear scaling. It is thus to be expected, that wavelet 
based techniques will catalyze progress in many fields of science and engineering, where 
such problems exist. An detailed tutorial style book describing how to use wavelets for 
the solution of partial differential equations will soon be published by the authors. 
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